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ABSTRACT 

In this paper we systematically study the spectrum and structure of incompressible MHD turbu- 
lence by means of high resolution direct numerical simulations. We considered both balanced and 
imbalanced (cross-helical) cases and simulated sub-Alfvenic as well as trans- Alfvenic turbulence. This 
paper extends numerics preliminarily reported in Beresnyak & Lazarian 2008. We confirm that driven 
imbalanced turbulence has a stationary state even for high degrees of imbalance. Our major finding 
is that the structure of the dominant and subdominant Alfvenic components are notably different. 
Using the most robust observed quantities, such as the energy ratio, we believe we can reject several 
existing models of strong imbalanced turbulence. 

Subject headings: MHD - turbulence ~ ISM: kinematics and dynamics 
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1. INTRODUCTION 

Astrophysical fluids are known to be turbulent from 
large scales of the order of hundreds of kiloparsecs, as 
in the case of galaxy clusters (Vogt & Ensslin 2005), to 
scales as small as hundreds of kilometers. Armstrong 
et al. (1995) reported a wide range of electron den- 
sity fluctuations probed by scintillation and scattering 
techniques. New statistical techniques using Doppler- 
shifted lines have enabled studies of velocity fluctuations 
on scales larger than a fraction of a parsec (see Lazarian 
(2009) and refs. therein). Recent years have been marked 
by new understanding of the key role that turbulence 
plays in a number of astrophysical processes (Cho, Lazar- 
ian & Vishniac 2003, Elmegreen & Scab 2004). Most no- 
tably, turbulence has drastically changed the paradigms 
of interstellar medium and molecular cloud evolution 
(Stone et al 1998, Ostriker et al. 2001, Vasquez-Semadeni 
et al. 2007, see also review by McKee & Ostriker 2007). 
Small scale turbulence has been probed by a variety of 
approaches such as gyrokinetics, Hall MHD and electron 
MHD (Howes et al 2006, Schekochihin et al 2007, Galtier 
et al 2003, Cho & Lazarian 2004). This progress calls for 
better understanding of the fundamentals of turbulence. 
One reason for doing this is to understand to what extent 
the turbulence simulated with diffusive computer codes 
resembles actual turbulence in astrophysical fluids with 
low diffusivity. 

We start with a few general remarks on the relation 
of fluid dynamics and description of astrophysical fluids. 
The continuous fluid description has been successful de- 
scribing a wide range of physical phenomena. In space, 
due to the presence of ionizing radiation and cosmic rays 
(CRs) the medium is almost always ionized to some de- 
gree, i.e. is a plasma. Although the dynamics of plasma 
is complicated on small scales, it can be considered well- 
conducting continuous fluid on large scales and, there- 
fore, Magnetohydrodynamics or MHD is applicable. 

One of the most interesting phenomena in fluids is 
turbulence, a seemingly random stochastic flow which 
appears spontaneously as long as the microscopic dissi- 
pation coefficients such as viscosity or magnetic diffu- 

Electronic address: [andrey, lazarian@astro.wisc.edu] 



sivity are small (which correspond to large Reynolds or 
magnetic Reynolds numbers). Turbulence increases dis- 
sipation due to so-called turbulent cascade, a nonlinear 
transfer of energy to smaller scales. As astrophysical flu- 
ids are turbulent, this affects dynamics and is manifested 
in a variety of situations such as reconnection, momen- 
tum transfer in accretion disks, etc. 

The basic theoretical study of the complicated non- 
linear dynamics of turbulence has been concentrated on 
incompressible turbulence with large Reynolds number 
(Kolmogorov, 1941). While dissipation in plasma might 
be much more complex than dissipation in molecular 
gases^, in the asymptotic large Reynolds number flows 
the "inertial range" fluctuations do not feel the peculiar- 
ities of the dissipation. On the other hand, it is often 
necessary to account for compressibility of the fluid, as 
it could be significant. An example of which is the inter- 
stellar medium which has has sonic Mach numbers in a 
range of 0.1 — 10. In this case the dynamics of Alfven and 
slow-mode perturbations on small scales can be consid- 
ered incompressible (Lithwick & Goldreich 2001), while 
fast mode has the dynamics of its own (Cho & Lazar- 
ian 2002, 2003). On the other hand, even at large scales 
where the compressibility is significant, one can find the 
features of incompressible dynamics (see, e.g., Beresnyak, 
Lazarian & Cho 2005). In other words, the understand- 
ing of incompressible MHD turbulence, which is often 
called Alfvenic turbulence due to the dominant role the 
Alfven perturbations play in the cascade, is of utmost 
importance. 

Ideal incompressible MHD equations can be written in 
the following simple form 



5tw=^ -I- S{-w^ ■ \/)w^ 



0, 



where 5 is a solenoidal projection operator and Elsasser 
variables are defined in terms of velocity v and magnetic 
field in velocity units b = B/(47r/9)^/^ as w+ = v -I- b 
and = V — b. Although fairly idealized, the problem 
of incompressible MHD turbulence with large Reynolds 
number has been difficult. First attempts to treat it was 

^ The dissipation in collisionless plasmas has been an area of 
active research recently, see, e.g. Howes et al 2008 and ref. therein. 
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an IK model by Iroshnikov (1963) and Kraichnan (1965). 
They found that the mean magnetic field provided by 
large scales will, very much unlike hydrodynamics, cru- 
cially define the dynamics on small scales. Although this 
first model was isotropic, it was eventually realized that 
the dynamics will result in anisotropy (Montgomery & 
Turner 1981, Shebalin ct al 1983, Higdon 1984). This 
resulted in Goldreich & Sridhar 1995 (henceforth GS95) 
model which postulated so-called "critical balance", i.e., 
maximum anisotropy which is allowed under strong in- 
teraction. At the same time, the concept of the dominant 
perpendicular cascade, that was used in GS95, has been 
validated in an analytical pcrturbativc theory of so-called 
weak Alfvcnic turbulence (Gahicr ot al. 2000, 2002)^. 

While hydrodynamic turbulence have only one energy 
cascade, the incompressible MHD turbulence has two, 
due to the exact conservation of the Elsasser (oppositely 
going wave packets') "energies". This can be also for- 
mulated as the conservation of total energy and cross- 
helicity^. The situation of zero total cross- helicity has 
been called "balanced" turbulence as the amount of op- 
positely moving wavepackets balance each other, the al- 
ternative being "imbalanced" turbulence. Most of the 
above studies concentrated on the balanced case, and, 
without exception, the GS95 model, which is the strong 
cascading model with critical balance, can only be kept 
self-consistent assuming balanced case. The real life tur- 
bulence, however, is almost always imbalanced, such as 
in situations when one has a strong localized source of 
perturbations (the Sun for solar wind or central engine 
for AGN jets), but also due to inhomogeneity of energy 
sources for turbulence (supernovas and stellar winds in 
the ISM) and the tendency of the decaying turbulence to 
become increasingly more imbalanced with time. More- 
over, the purely balanced MHD Alfvenic turbulence can 
not be understood as it is, without understanding of 
the more general imbalanced case. This is due to the 
fact that turbulence is a stochastic phenomena with all 
quantities fluctuating and every piece of turbulence at 
any given time can have imbalance in it. In this respect, 
while the mean-field Kolmogorov model can be expanded 
to include intermittency, the mean field GS95 model can 
not. 

Imbalanced turbulence, or "turbulence with non-zero 
cross-helicity" has been discussed long ago by a number 
of authors (Dobrovolny et al. 1980, Matthaeus & Mont- 
gomery 1980, Grappin et al. 1983, Pouquet et al. 1988, 
Biskamp 2003 and refs. therein). These papers testi- 
fied that the non-zero cross-helicity modifies the turbu- 
lence. Although these studies correctly reproduced sep- 
arate cascades for energy and cross-helicity, they were 
based on then-popular models of MHD turbulence and 

^ Although a successful analytical theory, the weak Alfvenic tur- 
bulence can rarely be applied to real-life MHD turbulence as it 
would require weak, relatively isotropic fluctuations on outer scale 
and a strong mean field. Even when these conditions are satis- 
fied, the strength of the interaction, which can be estimated as 
SviL/vaI, where L/l is the anisotropy of the perturbation, will 
increase down the cascade and break the applicability of the the- 
ory. In real life ISM turbulence the perturbations of the magnetic 
field are of the order of the mean field which makes the turbulence 
strong from outer scale inwards. For phenomcnological treatment 
of weak Alfvenic turbulence, see also Lazarian & Vishniac 1999. 

3 The latter, / v B X is a, quantity conserved in the absence 
of dissipation. 



later it became evident that these models arc problem- 
atic. For example, the closure theory of isotropic turbu- 
lence (Pouquet, Frisch & Leorat 1976) which reproduced 
IK model can be rejected by both theory and numer- 
ics'^. Another class of models were based on so-called 
two-dimensional MHD turbulence that, as we now know, 
is unable to reproduce important properties of the real 
three-dimensional turbulence, such as critical balance. 

Recently several models for the strong imbalanced tur- 
bulence have been proposed (e.g., Lithwick, Goldreich 
& Sridhar 2007, Beresnyak & Lazarian 2008, Chandran 
2008). As long as the full self-contained analytical model 
for strong turbulence continues to elude discovery, direct 
numerical simulations (DNS) will be an inspiration and 
guidance to theorists. While the Reynolds numbers in 
those simulations are fairly modest (800-4000 are the best 
to-dato), some of the robust quantities can be measured 
and used as a guidance to reject theories. While in this 
paper we are fully aware of these limitations, we present 
the most robust statistical measures, such as total ener- 
gies, dissipation rates and second-order structure func- 
tion (or its equivalent, the spectrum) and the anisotropy 
derived from it. We leave the study of higher order statis- 
tics (and intermittency) to the future studies. Our paper 
is written on the premise that one cannot fully confirm 
a model using three-dimensional DNS due to their fairly 
modest resolution, but from these DNS one can collect 
enough numerical evidence to reject a model. 

This paper expands parameter space of the prelimi- 
nary numerical results reported in Beresnyak & Lazar- 
ian (2008a). A short introduction into theories of im- 
balanced turbulence is in §2. Numerical code, the sim- 
ulation setup, and the establishment of the stationary 
state are explained in §3. The discussion of the structure 
function calculated parallel to the magnetic field (par- 
allel SF), a quantity which is a key for understanding 
statistical anisotropy of MHD turbulence is given in §4. 
In this section we explain various methods to define the 
local guiding field and the influence to the measurement 
of the parallel SF. Scale-dependent anisotropies derived 
from structure functions as well as power spectra are de- 
scribed in §5. Polarization alignment briefly discussed 
in §6. Final comparison with models as well as obser- 
vational data are in §7. We summarize our findings in 
§8. 

2. THEORETICAL CONSIDERATIONS 

The original GS95 model was based on the renormal- 
ization rule for Alfven wave's frequency called "criti- 
cal balance" . The necessity of such renormalization can 
be seen from a rigorous theory of weak Alfvenic turbu- 
lence (Galtier et al 2000, 2002) that predicts so-called 
"perpendicular cascade", i.e. the result that nonlin- 
ear interaction of Alfvenic waves conserve these wave's 
frequencies and only transverse structure of the wave 
packet is affected. As perpendicular cascade proceeds 
to small scales, the applicability of weak interaction 
breaks down, and Alfvenic turbulence becomes strong. 
In this situation GS95 assumed that the frequency of 

* The artificial term for "relaxation of triple correlations" , that 
was necessary to uphold local isotropy in this model, happen to 
be larger than real physical nonlinear interaction. Also numerics 
show that MHD turbulence is locally anisotropic (Cho & Vishniac 
2000, Maron & Goldreich 2001, Cho, Lazarian & Vishniac 2002) 
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the wavepacket can not be smaller than the inverse life- 
time of the wavepacket, estimated from nonlinear inter- 
action. In their closure model GS95 have an explicit 
ad-hoc term that allows for the increase of the wave 
frequency. The scale-dependency of this term is based 
on the assumption of turbulence locality (i.e. there is 
one characteristic amplitude of perturbation pertaining 
to each scale). In the imbalanced case, however, we have 
two such characteristic amplitudes and the choice for fre- 
quency renormalization becomes unclear (GS95)^. Any 
theory of strong imbalanced turbulence, which is quali- 
fied for serious consideration, must deal with this prob- 
lem. Let us first demonstrate that a direct generaliza- 
tion of GS95 for imbalanced case does not work, namely 
if we assume that the frequency renormalization for one 
wavepacket is determined by the amplitude of the oppo- 
sitely moving wavepacket. Indeed, in this situation the 
wave with small amplitude (say, w~) may only weakly 
perturb large amplitude wave w+ and the frequency of 
cascaded will conserve. On the other hand, may 
strongly perturb w~ and i/;~'s frequency will be deter- 
mined as /l^. This creates an inconsistency of the 
local cascade where both wavepackets must have com- 
parable parallel and perpendicular scales. In order to 
deal with this fundamental inconsistency, a new physical 
assumption must be adopted. Due to this fact, in the 
paper we mostly discuss three models of strong imbal- 
anced turbulence, Lithwick, Goldreich & Sridhar (2007) 
(LGS07), Beresnyak & Lazarian (2008a) (BL08a), Chan- 
dran (2008) (C08), that clearly state: a) the problem 
above; b) the new physical assumptions being adopted; 
c) an internally consistent physical model that follows 
from these assumptions and d) the full predictions of 
the turbulence spectra and anisotropy. In the discussion 
section of this paper we also mention other models that 
claim to have predictions on the strong imbalanced tur- 
bulence. We test these incomplete models by comparison 
with our numerics whenever possible. 

2.1. LGS07 and COS models 

LGS07 resolve the inconsistency explained above by as- 
suming that the strong wave is also cascaded strongly 
and its frequency is simply equal to the frequency of the 
weak wave, i.e. the critical balance for strong wave uses 
the amplitude of the strong wave itself (w+A = va^)- In 
other words, the anisotropics of the waves are identical. 
The formulae for energy cascading are strong cascading 
formulae, i.e. 

A 

This lead to the prediction /w^ = e+/e~. However, 
this prediction, together with assumption of the locality 

^ We assume that imbalanced turbulence is "strong" as long as 
the applicability of weak Alfvenic turbulence breaks down. This 
requires that at least one component is perturbed strongly. In the 
imbalanced turbulence the amplitude of the dominant component 
is larger, so that in the transition to strong regime the applicability 
of weak cascading of the subdominant component breaks down 
first. 

^ Throughout this paper we assume that 10+ is the larger- 
amplitude wave. This choice, however, is purely arbitrary and 
corresponds to the choice of positive versus negative total cross- 
helicity. 
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Fig. 1. — Upper: a w+ wavepacket, produced by cascading by 
w~ wavepacket is aligned with respect to w~ wavepacket, but 
misaligned with respect to the local mean field on scale Ai, by the 
angle d. Lower: the longitudinal scale A of the wavepackets, as 
a function of their transverse scale, A; A+, A~, Ai, A2 are the 
notations used in this paper. Modified from BLOSa. 

of cascading, lead to a contradiction on viscous scale^ 
where the nonlinear cascading rates must smoothly tran- 
sit into viscous dissipation rates. This require w"*" — w~ 
on the dissipation scale, so called "pinning" (Lithwick & 
Goldreich 2003). 

COS is a complicated quantitative theory of advection- 
diffusion cascading that have several rules to determine 
the diffusion of the waves in the parallel direction, which 
are analogous to the frequency renormalization of GS95. 
In effect, in the case of strong turbulence, the C08 rules 
lead to equal or very close anisotropics for both waves. 
Unlike LGS07, however, COS does not have a strong cas- 
cading for both waves, but (again, in effect) it has: 



X 



where the coefficients depend on the spectral slopes 
of w^. Although the theory of COS do not explicitly as- 
sume local cascading, in effect it produces such a locality 
as long as the ratio of dissipation rates (energy fluxes) 
e+/e~ is not very large (the critical value is around two), 
therefore COS also requires pinning on the viscous scale. 

2.2. BL08 model 

BLOSa relaxes the assumption of local cascading for 
the strong component w"*", while saying the w~ is cas- 
caded in a GS95-like way. In BLOSa picture the waves 
have different anisotropics (see Fig. [T]) and the w+ wave 
actually have smaller anisotropy than w~, which is oppo- 
site to what a naive application of critical balance would 
predict. The anisotropics of the waves are determined by 



w+(Ai)A {Xi)^vaXi, 
w+{X2)A+{X*)^vaXi, 



(1) 
(2) 



where A* = VA1A2, and the energy cascading is de- 
termined by weak cascading of the dominant wave and 
strong cascading of the subdominant wave: 

^ LGS07 does not discuss transition to viscous scales. 
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e+ = 



{w+{X2)fw-{Xi) w-(Ai)A-(Ai) 



^_ _ {w-{X^)fw+{M) 
Ai 



•/(A1/A2), (3) 



(4) 



One of the interesting properties of BL08a model is 
that, unUke LGS07 and COS, it does not produce self- 
similar (power-law) solutions when turbulence is driven 
with the same anisotropy for w+ and w~ on the outer 
scale. BLOSa, however, claim that, on sufficiently small 
scales, the initial non-power-law solution will transit into 
asymptotic power law solution that has A^/Aj = e"'"/e~ 
and A2/A1 = (e+/e-)3/2. The range of scales for the 
transition region was not specified by BLOSa, but it was 
assumed that larger imbalance will require larger transi- 
tion region. 

3. NUMERICAL SETUP 

Incompressible MHD equations with dissipation and 
driving are, 



5tw± + ^(w^F • V)w± = -z/„(-V2)"w± + f±, (1) 

where n is an order of hyperdiffusion, and is the 

driving force, whose rms values in arbitrary units are pre- 
sented in Table 1. These equations were solved by the 
code which is similar to one in Cho & Vishniac (2001). 
The differences include introduction of Elsasser driving, 
the handling of arbitrary physical sizes of the box, re- 
gardless of numerical resolution {\.q. elongated boxes for 
sub-Alfvenic turbulence) and significant improvements in 
numerical efficiency. Om code is pseudospectral, i.e., it 
solves ODE in time with finite difference for each spaeial 
Fourier harmonic, the nonlinear term being calculated 
in real space, while solenoidal projection and dissipa- 
tion terms were applied in Fourier space. Pseudospectral 
methods in fluid dynamics has been known since 80-s 
(Canuto et al, 19S8). The strong point of pseudospec- 
tral codes is that they allow precise control over dissipa- 
tion and exact incompressibility. Not only does it relieve 
worries with respect to grid effects and numerical dissipa- 
tion, but it also makes possible the use of hyperviscosity 
and hyperdiffusivity which extends useful inertial inter- 
val although at the expense of increased (compared to 
physical viscosity) bottleneck effect and a different form 
of the spectral slopes (for more discussion on spectral 
slopes, see Beresnyak & Lazarian 2008b (BL08b)). The 
simple version of our pseudospectral code uses periodic 
boundary conditions. This necessitates a discussion on 
whether this introduces artificial effects. In the subse- 
quent sections we show that a) our numerical boxes have 
enough parallel size to allow for eddies of the largest par- 
allel size, dictated by dynamics, to exist; b) at any given 
time our box contains large number of independent tur- 
bulent realizations {'^ 40) c) the dynamical time of the 
eddy is several times smaller than it takes the eddy to 
cross the box boundaries. 

3.1. Choice of physical dimensions, numerical 
resolution and driving 

The sub-Alfvenic turbulence, where the perturbation 
strengths are smaller that va is either weak or strong 




Fig. 2. — Two-dimensional energy spectrum in sub-Alfvenic and 
trans- Alfvenic case. For sub-Alfvenic case the abscissa kx is scaled 
with \/Ma factor. The dashed line indicates GS95 fcy ~ k^J^^ 
anisotropy. 

but anisotropic. The critical anisotropy is determined by 

the breakdown of the applicability of weak MHD turbu- 
lence which happens when k\^]^VAlk_\_5v ~ 1. In most of 
our simulations we drive turbulence on outer scale with 
the same anisotropy for both wave species, so this break- 
down is determined by 5v of the strong wave. We drive 
turbulence on outer scale in such a manner that the 
strong interaction establishes on the outer scale. Also, 
we provided driving for k = 2.. 3. 5 which means that 
the maximum eddy size is several times smaller that the 
box. This is to ensure that the first turbulent scales 
fc « 4 have more than enough space in parallel direc- 
tion in case that we did not estimate the transition into 
strong interaction regime correctly and the parallel scale 
of the cascaded eddies is longer than we expected. The 
results from §5, however, suggest that our choice was 
correct, with the largest coherent eddy size being around 
1/4 of the box size in both parallel and perpendicular di- 
rections. We used fully predetermined stochastic driving 
in both Elsasser variables * with a certain amplitude of 
the force (/^, see Table 1), so the energy input wasn't 
strictly controlled by the forcing, but rather was calcu- 
lated during the simulation^. In addition, we developed 
a driving which ensures constant energy input for both 
components, these tests confirm properties of the imbal- 
anced turbulence that were obtained with fully stochastic 
driving. We used the latter for most of our simulations. 

In studying sub-Alfvenic turbulence we adopted the 
approach to increase va by increasing Bq and increase 
the parallel physical size of the box L by the same fac- 
tor I/Ma without changing the equilibrium value of 5v, 
so that strong interaction timescale \/5v stays constant 

Elsasser driving is a preferred way to study intcrtial range of 
subAlfvenic turbulence as it simulates the supply of Elsasser en- 
ergies from larger eddies of a realistic turbulence. It is important 
to remember that kinetic and magnetic energies are not separately 
conserved by MHD equations. So when one has a pure velocity 
driving in a simulation with mean field (as in Cho & Vishniax; 
2001), he will generate approximately as much magnetic pertur- 
bations due to the Alfven effect, the result being two Alfven or 
pseudo-Alfven waves propagating in opposite directions. These 
waves, however, would have an artificial correlation (reflected by 
the fact that at t = b = 0). In order to use all degrees of free- 
dom and have better stochasticity one has to drive w+ and w~ 
independently. The mechanisms by which the outer scales of a re- 
alistic, say, ISM turbulence are driven are briefly discussed in §7.3 
and §7.4. 

^ Some simulations of hydrodynamic turbulence used negative 
viscosity on large scales to drive turbulence. In MHD this some- 
what unphysical approach docs not work because in this case even 
in the balanced simulations imbalance spontaneously occur and 
continue to increase without limit. 
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and similarly the eddy transverse time K/va also stays 
constant. Alternatively, one can keep -Bo constant, but 
decrease 5v^ but in this case the timescales of sub- and 
trans- Alfvenic turbulence will be different. Also it is 
harder to tune the equilibrium 5v, rather than Bq and 
L. 

Note, that one can naively assume that due to GS95 
anisotropy one needs lower numerical resolution in the 
parallel direction, approximately by the ratio of the 
anisotropics on the driving scale and on the dissipation 

scale, which is {k±diss / ki_drivY^^ in the GS95 model, and 
can be a number between 2 and 4 in a high resolution 
MHD simulation. For instance. Bigot, Galtier & Politano 
2008 used 512x512x64 numerical resolution. On the sec- 
ond thought this approach is not evident, since the high- 
est values of in the global reference frame will be de- 
termined by field wandering on the outer scale. In other 
words the anisotropy in the global frame will be approxi- 
mately scale-independent and the ratios of k±diss/k±driv 



and fcii 



[diss / k^^driv 



will be almost equal, which necessitates 
the use of NxNxN numerical resolution, i.e. cubes, for 
both elongated {AIa < 1) and cubic {Ma ~ 1) physical 
boxes. 

We confirmed this by plotting the parallel and perpen- 
dicular spectra in the global frame and saw that the par- 
allel spectrum protrude to almost as far as k±rnax/MA- 
Fig. [2] shows how energy is distributed on the two- 
dimensional fc||,fcj^ plane (global reference frame). We 
see that while most of the energy is in GS95 cone, there 
is also plenty of energy outside of it, especially in the 
upper right corner which correspond to maximum space 
frequencies in both parallel and perpendicular direction. 
If one decide to significantly cut numerical resolution in 
parallel direction he/she would incorrectly describe the 
dynamics on small scales. In only one of our simula- 
tions, Al (see Table 1) we were able to cut parallel res- 
olution by a moderate factor of 1.5 without sacrificing 
small parallel scales, due to the relative lack of energy in 
parallel direction in this particular balanced sub-Alfvenic 
case. In all other simulations such a reduction was not 
possible because most of the k-space was filled with en- 
ergy. We note that Miiller & Grappin (2005), by using 
1024x1024x256 resolution in their balanced sub-Alfvenic 
simulations have reduced parallel resolution by a factor 
of 4, which is, most likely, too large. 

For all simulations A1-A8 we used hyperviscosity and 
hyperdiffusivity of 6th order (fc^). This choice was ne- 
cessitated by the nature of imbalanced turbulence which 
has shorter inertial range for dominant wave due to fairly 
large cascading timescale of this wave (see §2). With 
currently available numerical resolutions one cannot see 
an inertial interval of the strong wave in a simulation 
with large imbalance and real (fc^) diffusivity. Unfortu- 
nately, due to the bottleneck effect, hyperdiffusion have 
affected spectral slopes, although the effect on anisotropy 
was much less. We refer to BLOSb for a comparison 
of turbulent simulations with normal and hyperviscos- 
ity. Due to hyperviscosity the dissipation scale was fairly 
small, the dissipation cutoff was around k = 200 (with 
Nyquist frequency of 384) for balanced simulations and 
about the same for weak component in imbalanced simu- 
lations. The strong component for the most imbalanced 
simulations A7 and A8 had a cutoff around k — 100 
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Fig. 3. — The energy evolution in time for both Alfvenic modes 
for e+/e- = 12. 

(Fig. [5]). Due to hyperviscosity we can not uniquely de- 
fine a Reynolds number of our simulations, however vis- 
cous simulations with Re = Rem ~ 6000 could provide 
turbulence inertial ranges that are similar to ours. 

Once we have chosen the geometry of our simulation 
and figured out the extend of the perturbations on the 
spectral plane, the choice of timestep becomes evident. 
On one hand, for the dissipation term we use integra- 
tion technique (Cho & Vishniac 2000, see also Maron 
& Goldreich 2001), and since we don't worry too much 
about the precision of the dissipation term, it doesn't 
limit the timestep. On the other hand, the general non- 
linear term, containing both Bq and 6v can be seen as 
the sum of linear advection term with velocity va and 
nonlinear advection with Sv, 5b, etc. In turbulence, that 
is driven to be strong on the outer scale, these terms 
will be of the same order if we refer to the outer scale, 



i.e. the terms will be VA5vk±driv and Sv6vk 



\driv • 



On 



the dissipation scale these terms will be determined by 
VA^vkji^diss and 5v5vk\\diss^ which are, by the argument 
in the above paragraph, again on the same order. So we 
can just use linear advection behavior to estimate the 
timestep. This behavior in k-space is, essentially, a rota- 
tion of the phase of the wave, in a manner of exp{ik\\ VAt) . 
In order to reproduce this rotation numerically we need 
k\\maxVASt to bc Smaller than unity, such as around 0.1, 
so that the code stays stable, since we don't need good 
precision beyond the dissipation scale where there is no 
energy. 

The average dissipation rates reported in Table 1 
were calculated using a sum of the work done to the 
Elsasser fields, i.e. we summed (w='= -I- f^dt) ■ i^dt at 
every timestep. As our code (its nonlinear part) was en- 
ergy conserving, we assume that the same amount of 
energy was, on average, lost to the dissipation term. 
We also confirmed these values of by using third 
order Chandrasekhar-Politano-Pouquet structure func- 
tions (see, e.g., Biskamp 2003), which quantify nonlinear 
energy transfer. 

3.2. Establishment of the stationary state 

One of the goals of this paper is to demonstrate that 
a stationary state exists for imbalanced turbulence with 
rather high degree of imbalance. Note, that the local 
model of weak Alfvenic turbulence work for imbalances 
of no more than e+/e~ = 2 (Galtier 2000, Lithwick & 
Goldreich 2003), and the model of strong imbalance tur- 
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768^ 


10:1:1 
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0.16 


0.08 


154 


33 


12.0 


6 


2 


7.4 


145 ± 10 


A6 


7683 


1:1:1 


1 


0.19 


0.095 


40 


33 


12.0 


5 


2 


5.4 


90± 10 


A7 


768^ 


10:1:1 


10 


0.06 


0.02 


204 


63 


21.0 


8 


3 


16 


1150 ± 100 


A8 


768^ 


1:1:1 


1 


0.075 


0.025 


160 


63 


21.0 


8 


3 


12 


1100 ±100 



TABLE 1 

Simulations of strong sub-Alfvenic (Bq = 10) and trans-Alfvenic (Bq = 1) turbulent flows. A1 and A2 are balanced 
simulations for comparison with the rest. a3 and a4 are slightly imbalanced, a5 and a6 are strongly imbalanced and a7 
and a8 are very strongly imbalanced. aio is the duration of prior low resolution run in code units, at is the duration of 
the high resolution run, interval atl in the end of high-resolution run was used for data analysis. 




k 

Fig. 4. — Spectrum relaxation towards stationary state. The 
resolution of the experiment was increased from 128^ to 256'^. Only 
the spectrum of the dominant wave is shown. The lines correspond 
to t = 15,30,60, 150. 

bulence of COS also require similar limitation. 

The highest imbalance we attempted in our simulation 
of £"*"/€" = 16 was essentially limited by the long times of 
establishment of the stationary state. Note that accord- 
ing to BLOSa the dominant wave is cascaded weakly and 
its cascading times could be very large. Fig. [3] shows the 
total energy evolution for both modes for the e^/e^ = 16 
case^". The full relaxation towards stationary state re- 
quired around 300 Alfven times or 50 crossing times. 

As high imbalanced simulations proved to be so com- 
putationally expensive, we made a second experiment, 
which was to take the initial state that was already sta- 
tionary and to increase the numerical resolution, which 
allow the spectrum to extend to larger wavenumbers. 
Note, that our forcing, although stochastic, was prede- 
termined for each particular simulation and did not de- 
pend on numerical resolution. Now the question was how 
fast the spectra will relax to their stationary states. It 
turned out that the spectrum of the sub-dominant wave 
relaxed almost instantly, in a one dynamical time, which 
is consistent with BLOSa, while for the dominant wave 
the relaxation time was long. Note, the the dynamic 
(kinetic) timescale l/v for this region of /c-space for the 
strong wave was rather small, around 0.3. The relax- 
ation process is shown on Fig. [H It took 5t k, 60 to 

Time was measured in Alfvenic units, but the size of the box 
was 27r, thus 27r was the time for an eddy to cross the box. The 
time for the largest turbulent eddy to cross itself, and also the 
largest eddy dynamical time (L/v) was around unity, because the 
size of the largest eddy was a fraction of around 0.2 or 0.3 of the 
cube size (see Figs. [7]and|8ll. 
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Fig. 5. — The ratio between parallel SFs calculated with different 
definition of the local field. Upper plot is for balanced simulation, 
while two lower plots are for strongly imbalanced simulation. The 
reference SF is "minimum" one described in the text. The dotted 
line is "model-dependent" SF, dashed line is for "following a field 
line" method and the dotted-dashed line is for constant global mean 
field. 

get reasonably close to the stationary state. We consid- 
ered this experiment a success and used the technique of 
increasing resolution to save computational time. 

We also studied long-term evolution of nearly-balanced 
case when we allowed the low-resolution version of A3 to 
evolve for 500 time units. We didn't notice any long- 
term trends either in three-dimensional spectrum or in 
any other quantities during this run. 

4. PARALLEL STRUCTURE FUNCTION 

Unlike perpendicular structure function which is 
largely insensitive to the direction of the local field, the 
definition of the local field strongly affect the parallel 
structure function, which, in turn, determines the shape 
of the turbulent eddy. Since the latter is the major object 
of study in this paper we feel that the proper explanation 
of this point is due. 

The simplest way to define parallel structure function, 
S'J^II , is to take samples along the global mean field. This 
definition is, however, fairly bad, as it does not take into 
account field wandering. We expect SF\\ , defined in such 
way, along with perpendicular structure function to re- 
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0.01 0,10 

perpendicular scale 

Fig. 6. — This plot shows values of Aavr at which the mini- 
mum of the parallel structure function is reached. Triangles show 
perpendicular scales Aavr at which the minimum of SF^{w~,A) 

is reached, while squares show perpendicular averaging scales at 
which the minimum of SF^{w^,A) is reached. Solid and dotted 

lines indicate and w~ eddies' anisotropy which are defined in 
§5 and presented on Fig. 1111 

fleet the anisotropy in the global frame, which, by the ef- 
fects of field wandering, as we argued in §2, will be similar 
to the outer scale anisotropy. Therefore, such definition 
will effectively erase scale- dependent anisotropy which is 
the property of GS95-type models. 

Another way is to define local magnetic field by averag- 
ing over some scale A. In this way the parallel structure 
function become a function of two scales, such as 

SF^{w^,A, A) = {{w^{r - Ah^/b^) - w^{r))^),, 

where is the magnetic field averaged over scale 
A. We use Gaussian averaging defined as — 
l/XV2^ /b(r-R)exp(-i?V2A2)dR. In order to re- 
duce such a SF to a function of only A one can introduce 
a dependence between A and A, and plug in the A = /(A) 
in the above equation. This definition of SF^^ will be a 
model- dependent though. 

5F|2_,,,(z«±,A) = 5F|f(w;±,A,/(A)) 

For the balanced turbulence the anisotropy was mea- 
sured to be close to the one predicted by GS95, i.e. 
A ~ A^/^. One can, therefore, introduce a reasonable 
model-dependent SF^^ as taking f{x) — const ■ x^l"^, 
where the constant depend on the outer scale of the sim- 
ulation. As we show below, this definition is almost per- 
fect for balanced turbulence, but the question is whether 
it does equally well for the imbalanced case. 

Let us consider some model-independent ways to deter- 
mine SF\\ . Apparently, the first definition using global 
field is model-independent, but fairly bad, it correspond 
to taking averaging A = oo. One can also take A = 0, 
i.e. always use local field without any averaging. An in- 
teresting model-independent method was used in Maron 
& Goldreich 2001, where two points were chosen to lie 
on the same magnetic field line. The distance A was also 
calculated along the line. 

When we look for anisotropy we normally want to ob- 
tain lower values of SF\^ . According to the eddy ansatz. 




Fig. 7. — Comparison of the SFs from trans- Alfvcnic (left) and 
sub-Alfvenic (right) balanced simulations. Note the difference in x 
axis between two plots which indicates that Al is approximately 
10 times more anisotropic. Contours indicate SF levels, solid line 
is a demonstration of GS95 A ~ A^/^ law 

outlined in §2, we receive Zower values of |w^(r— An(r))— 
w^(r)|, where n is a unit vector along the eddy. There- 
fore, the averaging of the field that provides minimum 
values of SF\^ approximates the direction of the eddy 
alignment better, provided that there is a connection be- 
tween the field direction and eddy alignment (if there is 
no such connection, there will be no dependence on the 
averaging scale A). 

So, another model independent way to define SF\^ will 
be 

SFl^^r. i^^ : A) = mm SF^ , A, A) . 

This definition not only provides us with the value of 
S'i^ll, but, giving A at which minimum is achieved, gives 
us a hint to how eddies are aligned with respect to the 
magnetic field. 

Fig. [5] shows a comparison between different meth- 
ods to calculate parallel SF. We plotted them relative to 
^F^min- balanced case the three methods - "min- 

imal" , "following the field line" and "model-dependent" 
work very well, while "global field" method doesn't work. 
The latter confirms that turbulent eddies are aligned 
with respect to local field, not the global field (Cho & 
Vishniac 2000). In the imbalanced case the situation is 
more complicated. For the weak component all three 
"good" methods work very well, while for the strong 
wave there is a systematic error for all methods, except 
"minimal" . This is due to the fact that in the imbal- 
anced turbulence the strong component (w~^) eddies are 
aligned with respect to much larger scales of the mag- 
netic field (see §2.2). Since most magnetic field pertur- 
bation is provided by the strong wave, it follows that the 
strong field eddies are aligned with their own field on 
a larger scale. This directly confirms the prediction of 
BL08a model. On the bottom panel of Fig. [5] the "field 
line" method gives values that are smaller than "mini- 
mal" method. This is due to the fact that in the field 
line method we measured the distance along magnetic 
field, and the physical distance was actually shorter, this 
allowed for smaller than SF^^-^ values on the outer scale 
where the difference between straight-line distance and 
along-the-field-line distance is significant. 

It turns out (Fig. [6]) that the averaging scale at which 
minimum of parallel structure is reached for weak wave 
approximately corresponded to its anisotropy, which is 
consistent with strong cascading hypothesis. But for 
strong wave this averaging scale is larger than the per- 
pendicular scale dictated by anisotropy. In other words. 
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Fig. 8. — Comparison of the SFs from dominant (left) and sub- 
dominant (right) Alfvcnic waves from = 16 imbalanced sim- 
ulations. The anisotropies of components are notably different. 
The solid lino is the GS95 A ~ A^/S law 

the eddies of the strong wave are ahgned with respect to 
the magnetic field which is averaged on larger scale than 
the eddy's own perpendicular scale. This is consistent 
with the BL08a model. 

5. SPECTRA AND ANISOTROPIES 

We calculated two-dimensional (depending on paral- 
lel and perpendicular distances) second order structure 
functions with respect to the local field using "model 
dependent" definition of the local field from previous 
section. Although this method slightly underestimates 
anisotropy, according to Fig. [S] it works fairly well. The 
structure functions were calculated using all available 
"stationary state" datacubes, i.e. were averaged over 
time. The contours of these SFs for balanced simulations 
Al and A2 are presented on Fig. [7] and for imbalanced 
simulation A7 on Fig. [51 Fig. [7] shows SFs for total 
energy i.e. it is summed over and . These fig- 
ures basically validate our assumptions from §3 regarding 
physical and computational dimensions of the box. We 
see that, according to expectations, trans- Alfvenic A2 is 
almost isotropic on outer scale but becomes progressively 
anisotropic towards small scales, while, as we expected, 
sub- Alfvenic Al has approximately 10:1 anisotropy on 
outer scale, and increases towards small scales. If we de- 
crease anisotropy of Al by a factor of 10 by rescaling 
X-axis we almost reproduce A2, the difference is mostly 
being on the outer scale (this difference is easier to see on 
Fig[Tni). Fig. [8] shows SFs for two separate components 
w+ and in the strongly imbalanced case of A7. The 
anisotropy on outer scale is approximately 10:1 for both 
components, which validates our choice of computational 
box. This anisotropy increases towards small scale, but 
in a different fashion for each component. We see the 
anisotropy of strong wave is almost 5 times smaller on 
dissipation scales. 

Fig.[9]shows so-called three-dimensional angle-summed 
spectra for both components in all simulations. These 
spectra are obtained by summation of spectra over solid 
angle for all wavevectors with the same magnitude k. It 
can be related to three-dimensional angle-averaged spec- 
tra by dividing by fc^. In the sub- Alfvenic cases Al, A3, 
A5 and A7, this spectrum is almost identical to the so- 
called perpendicular spectrum, which takes into account 
only structures perpendicular to the magnetic field and 
is the main target of prediction of GS95 model. As they 
are almost identical we did not have to plot perpendicu- 
lar spectrum separately. Another definition of spectrum 
which depend only on the magnitude of the wavevec- 
tor is the so-called one-dimensional spectrum (see, e.g. 
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Fig. 9. — Spectra for all data, gray shows mean-square fluctua- 
tions in time. 




0.01 0.10 0.01 0.10 

perpendicular distance 

Fig. 10. — Anisotropies for balanced simulations. 

Monin & Yaglom 1965). This spectrum is less sensitive 
to the bottleneck effect. We refer to the paper of Beres- 
nyak & Lazarian (2008b) (henceforth BL08b) for a more 
thorough comparison between one-dimensional and 3D 
spectra and discussion on bottleneck effect. 

In Fig. [21 the two bottom plots have relatively large 
variation (gray areas) on the end of the w'^ spectra. This 
is due to the fact that in A7 and A8 we barely reached 
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Fig. 11. — Anisotropies for imbalanced simulations. The mapping of A(A) is explained in 
and w~ increases with increasing imbalance. 



j5. The difference in anisotropy between 



stationary state (for more discussion, see Fig. |4]and §3.2) 
in the high-resolution run. 

Fig. [5] shows spectral slopes between -1.12 and -1.93 
with balanced simulations having slopes as flat as -1.37^^. 
The GS95 prediction is Kolmogorov's —5/3 ~ —1.67, 
while Boldyrev 2006 prediction is -1.5. The flat slopes 
observed in real data are most certainly due to rather 
strong bottleneck effect seen in simulations with hyper- 
viscosity. In the imbalanced case the predictions are fol- 
lowing: LGS07 predicts -1.67 slopes for all 8 cases; COS 
predicts -1.67 for balanced cases Al and A2, A5-A8 are 
outside of the applicability of his model, A3 and A4 must 
show very different slopes - approximately -1 for weak 
component and -3 for strong component, also COS pre- 
dicts pinning on dissipation scales i.e. spectra should 
converge on dissipation scale The ratios of the total 
energies (see Table 1) are predicted as following: COS - 
A5-A8 are outside of applicability of his model, A3 and 
A4 should have very large imbalance {w^)'^ j (w~Y of at 
least a 1000, while 4-6 is actually observed; LGS07 pre- 
dicted/observed - A3: 4/5.5, A4: 3/3.9, A5: 55/145, 
A6: 30/90, A7: 260/1150, AS: 144/1100. We see that in 
this respect deviations from LSG07 predictions are small 
for small imbalances but fairly large for large imbalances. 
BLOSa argues that if one drives turbulence with the same 
anisotropy on outer scale (as in these simulations) the 
anisotropies of the components will diverge towards small 
scales, and this solution will not be self-similar (and not 
power-law). However, BLOSa makes predictions regard- 
ing local slopes even in this case. This can be seen from 
(1) which is a classic critical balance between weak wave 
anisotropy and strong wave amplitude and (4) which is 
strong cascading of the weak wave. We do not expect 

The slopes for one-dimensional spectra are steeper, with bal- 
anced slope of -1.45 and the most imbalanced slopes -1.97 for strong 
wave and -1.22 for weak wave 

^"^ LGS07 does not discuss transition to viscous scales, but, as a 
model of local cascading, it must have pinning on viscous scale. The 
pinning, however, is impossible within the framework of LGS07, as 
the latter predicts to "'"/to" = 



relations between slopes based on (4) to hold, because 
it is strongly influenced by bottleneck effect (BLOSa also 
predicts -1.67 slopes for balanced case). However, there 
is some dependence between energy slope and anisotropy 
slope, similar to what BLOSa predicts. Namely, it follows 
from (1) that shallower anisotropy slope for means 
steeper spectral slope for which is observed. Also, 
from (4), steeper spectral slope for also means shal- 
lower spectral slope for w~ which is also observed. 

The anisotropy was measured in the following man- 
ner. First, the parallel and perpendicular second order 
structure functions were calculated, then we found equal 
values of parallel and perpendicular SFs and in this way 
the mapping or function between independent variables, 
parallel or perpendicular scales were created. This func- 
tion is plotted on Figs [TOl [TT] with shades of gray indi- 
cating RMS fluctuations in time. This definition of A(A) 
mapping can be understood from two-dimensional plot 
of SF, e.g. Fig. [71 when one follows a contour of SF 
and finds which parallel scale correspond to particular 
perpendicular scale. We see that for the imbalanced case 
anisotropy curves have different slopes and diverge from 
outer scale where they are equal (this is dictated by driv- 
ing) to smaller scales where they are different. 

We devoted §4 to the discussion of the measurements 
of the parallel structure function which was used in the 
above definition of the anisotropy curves. Although it 
might appear that each and every definition produce a 
different anisotropy curve, the major difference is be- 
tween global and local definition of the field direction, 
while all local methods ( "field line" , "model-dependent" 
and "minimal"), also dubbed "good" in §4, give very 
similar results. In fact, these is no perceivable quali- 
tative difference between anisotropy curves obtained by 
either local methods. This could be explaned by Fig. 4 
middle and bottom panels where the quantitative differ- 
ences between methods are small, but on the other hand, 
the dependence of SF\^ on scale is strong (~ Z||). Also, in 
the middle panel, the difference is mostly by a constant, 
which will only give a slight shift of the anisotropy plot. 
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Fig. 12. — Polarization angle alignment (dashed) and polar- 
ization intermittency (see definition in the text). Numbers are 
power-law slopes determined in the middle of the inertial interval. 

All in all, the claim that anisotopy curves will diverge 
by a factor of 3 to 4 in strongly imbalanced simulations 
stay true regardless of the "local" method used. We re- 
jected "global field" method as it does not reveal scale 
dependent anisotropy - a groundbase of GS95 model. It 
is worth noting that LGS07, COS and BL08 use GS95 
as a basis and smoothly transit to GS95 in the balanced 
limit. There is a wealth of theoretical arguments why 
the SFs has to be measured with respect to the local 
field (Cho & Vishniac 2000, Maron & Goldreich 2001, 
etc). We also would like to note that aside from driven 
simulations described in this paper we also observed a 
significant difference in vS^ and w~ anisotropics in a de- 
caying imbalanced simulations. 

COS and LGS07 both predict identical GS95 anisotropy 
for both modes, which is inconsistent with simulations. 
BLOSa predicts diverging anisotropy, most notably, with 
stronger wave having smaller anisotropy, which is consis- 
tent with simulations. The value of the differences, how- 
ever, do not reach the asymptotic value of j e~ which 
may be attributed to the short inertial range. 

6. POLARIZATION ALIGNMENT 

Aside from energy-type statistics (second order SFs 
and spectra) one can measure so-called alignment effects 
which are, in a sense deviations from the assumptions 
of independent randomness of fluctuations included in 
mean-field models. These were discussed in Boldyrev 
(2005) and numerically discovered for the first time in 
Beresnyak & Lazarian (2006). 

Fig 1121 shows two measures of alignment - the an- 
gle polarization alignment A A — (|sin6'|) (dashed line), 
where B is an angle between Elsasser variables perturba- 
tions i5w+ = (5v -|- (5b and (5w^ = (5v — (5b and "polariza- 
tion intermittency" PI = {\Sw^6w^ sin d\)/ {\Sw^6w~\) 
(solid line). 

For the more detailed numerical study of alignment 
effects we refer to Beresnyak & Lazarian (2006) and 
BLOSb. The potential effect of alignment on the energy 
cascade was discussed in Boldyrev (2005, 2006), although 
this effect has not yet been convincingly confirmed by nu- 
merics (see BLOSb for a discussion and refs). 



7. DISCUSSION 

7.1. Comparison with models 

We briefly mentioned the tentative nature of confirm- 
ing a model in §1, we also claimed the ability of numerics 
to reject some models on the basis of robust quantities. 
Although a theory can make a wide variety of predictions, 
only few of those can be effectively attacked by numerics. 
One of the quantities that is notoriously hard to measure 
in DNS is the spectral slope of turbulence. A difference 
between —3/2 slope and —5/3 can be masked by a va- 
riety of effects such as bottleneck effect, driving, and so 
on. In contrast, the quantities such as Kolmogorov con- 
stant are fairly easy to obtain and quickly converge with 
increasing resolution. In fact, modest resolutions such as 
12S^ give reasonably precise estimates of this constant. 
This is due to the fact that the total energy and the total 
dissipation rate are easy to measure, to get a statistical 
average, and also are free of uncertainties of interpreta- 
tion. What sort of models can be judged on the basis of 
these quantities? Such are the models of local cascading 
where the cascade rate depend only on the characteristic 
quantities of, say, wf^ on a particular scale / (and, pos- 
sibly, weakly depend on the spectral slope). In this case 
numerics only has to reproduce a one or two steps of such 
cascading to obtain reasonable dissipation rate based on 
a particular total energy. In this sense our testing does 
fairly well, as we mostly consider models of local cascad- 
ing (LGS07, COS, Perez & Boldyrev (2009), Podesta & 
Bhattacharjee (2009)), at the same time, with 768'^ nu- 
merical resolution we reproduce five to six binary steps 
in fc-space. 

Our numerical data strongly contradict to three mod- 
els of imbalanced turbulence, namely LGS07, COS, and 
Perez & Boldyrev (2009). In particular, two of these 
models COS, and Perez & Boldyrev (2009) show gross 
inconsistencies between observed and predicted energy 
ratios vs dissipation ratios. Indeed, COS must have a 
huge energy ratio (of around a 1000) in simulations with 
e+/e~ close to two (A3 and A4), while a modest ratios 
of 4 and 6 are observed. Perez & Boldyrev (2009) does 
extremely bad in cases with large imbalances. A7 and AS 
has an energy ratios of around a 1000, while predicted 
quantities are 16 and 12. LGS07 does a much better 
job on energy ratios, but still fails the A7 and AS (large 
imbalance) tests, see Fig. [T3l 

Furthermore, LGS07 and COS have predictions re- 
garding eddie anisotropics. Both of these models pre- 
dict equal anisotropics for and w~ while a different 
anisotropics are observed. Aside from these inconsisten- 
cies, we also note that COS predicts pinning at the dis- 
sipation scale, which is not observed. COS also predicts 
that the strong wave has to have steeper spectral slope 
than the weak wave. This corresponds to numerics quali- 
tatively, but not quantitatively, indeed, according to COS, 
A3 and A4 must have a huge slope difference of around 
2, while the real difference is around 0.12. 

Although it is harder to confirm a model rather than 
to reject a model by direct numerical simulations, we see 
that there is a qualitative agreement between BLOSa and 
numerics. Most of the features predicted by BLOSa are 
observed in simulations, namely a) the anisotropics of the 
waves are different and strong wave anisotropy is smaller; 
b) while the weak wave eddies are aligned with respect to 
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Fig. 13. — Total energy ratio versus dissipation rate ratio 
(see also Table 1). Diamonds: subAlvenic simulations, stars: 
transAlvenic simulations, solid line: LGS07 prediction, dashed line: 
Perez & Boldyrev (2009) prediction. There is no simple formula 
for COS model, but the energy ratio is expected to become very 
large when dissipation rate ratio approaches the critical value of 
around two. 

the local field on the same scale as the eddy, the strong 
wave eddies are aligned with respect to a larger-scale 
field, (Fig. (5]); c) the energy imbalance is higher than 
in the case when both waves are cascaded strongly (Ta- 
ble 1), which suggest that the strong wave is cascaded 
weakly; d) the dissipation scales for the weak and the 
strong waves are different, namely the inertial range for 
weak wave is longer (Fig. [9|) which is what predicted by 
BL08a; e) there is no "pinning" at the dissipation scale 
which suggest nonlocal cascading. 

We note that there is no quantitative agreement be- 
tween the difference in anisotropics in the "asymptotic 
power-law solutions" of BL08a and simulations (c.f. §2 
and Fig. [TT|) . This is probably due to the fact that 
asymptotic power-law solutions have not been estab- 
lished within our inertial range. 

Aside from three models with detailed strong imbal- 
anced theories that we described above, LGS07, BL08a 
and COB, recently two papers also made predictions re- 
garding strong imbalanced turbulence. Perez & Boldyrev 
(2009) predicted that (w+Zw")^ = which comes 

to much stronger contradiction with A7 and AS than 
LGS07 and as such should be discarded. We note paren- 
thetically that measurements of the stationary levels of 
energies for w"*" and w~ are the most robust and model 
independent of all other measurements. Another paper, 
by Podesta & Bhattacharjee (2009) is a modification of 
Perez & Boldyrev (2009) which tries to resolve a huge 
discrepancy of the latter. Unfortunately, their theory 
has an arbitrary parameter that can be tuned to change 
the prediction for (w^/w")^ and, therefore, it can not 
be tested numerically. We also note that two aforemen- 
tioned papers do not have clearly stated predictions for 
anisotropy, which also limits our ability to test them. 

7.2. Comparison with earlier simulations 

Maron & Goldreich (2001) and Cho et al. (2002) per- 
formed 3D numerical simulations of decaying imbalanced 
turbulence. In these simulations the perturbations ob- 
tained as a result of balanced driven MHD modeling were 
separated into oppositely moving flows of Elsasser vari- 
ables and one of the flows was arbitrary decreased in am- 
plitude. The authors observed the increase of the damp- 
ing time for the strong component. They also observed 
the increase of turbulence imbalance as the turbulence 



was evolving. Naturally, no stationary state was achiev- 
able for the imbalanced turbulence induced this way. 

7.3. Role of homogeneous turbulence driving 

Properties of imbalanced turbulence may depend on 
how it is driven. All three major models of imbal- 
anced turbulence that we discuss above (LGS07, BLOSa 
and COS) describe imbalanced turbulence driven homo- 
geneously through the volume. The properties of such 
turbulence may differ if the sources of driving are local- 
ized. To illustrate this fact, in the Appendix A we discuss 
a toy model of weak imbalanced turbulence driven inho- 
mogeneously at boundaries. Note that this approach is 
very different that the DNS approach of the rest of this 
paper. Being the model of weak turbulence, this toy 
model describes only a perpendicular cascade just like its 
precedessor Lithwick & Goldreich (2003), which consid- 
ered homogeneous case. However, the stationary states 
of these inhomogeneous and homogeneous cases substan- 
tially differ. We also expect to see substantial differ- 
ences for the homogeneous and inhomogeneous driving 
for the strong imbalanced turbulence, but since such a 
study with the use of the DNS will be complicated and 
computationally expensive, we defer this discussion to 
future publications. 

7.4. Effects of compressibility 

The simulations presented here are of incompressible 
MHD turbulence. Not only full compressible MHD equa- 
tions have more degrees of freedom (incompressible ap- 
proach exclude fast mode perturbations), but compress- 
ibility substantially changes the properties of turbulence 
when sonic Mach numbers are of around unity or higher. 
However, numerical studies in Cho & Lazarian (2002, 
2003) showed that coupling of Alfvenic and magnetosonic 
waves in strong turbulence is not as strong as was ex- 
pected, which can be due to the fast non-linear decay of 
Alfvenic eddies. 

In the imbalanced turbulence the strong wave is long- 
lived. Therefore, one can expect the imbalanced tur- 
bulence to be more affected by coupling of incompress- 
ible and compressible motions. Density inhomogeneities 
present in the compressible fluid can act as mirrors re- 
flecting waves and decreasing the degree of turbulence 
imbalance. Parametric instabilities (see, e.g., Del Zanna 
et al. 2001) can develop in the compressible fluid, de- 
creasing the imbalance. Thus, stationary states with 
high degree of imbalance may not be attainable in com- 
pressible fluids. Further research in this direction is nec- 
essary. 

8. SUMMARY 

In the paper above we have performed MHD numeri- 
cal simulations of the homogeneously driven imbalanced 
turbulence and have shown that: 

1. Stationary states exist for rather high degree of im- 
balance. 

2. For large imbalances, the ratio of the amplitudes dis- 
agrees with the predictions in LGS07. 

3. Rough correspondence of the expectations and mea- 
surements are obtained for the model in BLOSa, but more 
studies and testing is necessary. 

9. APPENDIX A: INHOMOGENEOUS WEAK IMBALANCED 
TURBULENCE - A TOY PROBLEM. 
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Fig. 14. — Spectrum of up (solid) and downgoing (dashed) waves, 
in the middle between the wave sources. 

Although we understand that the models of imbal- 
anced turbulence are in their making and even the case of 
homogeneous imbalance turbulence is still being investi- 
gated, we decided to add this appendix that probes into a 
different problem, namely, inhomogeneous Alfvenic tur- 
bulence. 

This is motivated by the fact that imbalanced turbu- 
lence often appears in inhomogeneous setting, i.e. when 
one has a strong localized source of waves, such as the 
Sun in solar wind turbulence. In this respect the the- 
ory of homogeneous imbalance turbulence should be well- 
applicable to small scales where the time scales are much 
smaller than outer timescales. However, on larger scales, 
homogeneity could be broken. In this section we probe 
the situation when it is broken. This can be used as 
a guidance to what extent the observations of the solar 
wind can be considered as restrictive measurements with 
respect to the theories of imbalanced turbulence. 

Previous studies of weak Alfvenic turbulence (e.g. 
Galtier et al 2000, Lithwick & Goldrcich 2003) considered 
spatially homogeneous case, which presume that and 
•w~ are driven homogeneously in space. This assump- 
tion is often violated in nature, specifically when we have 
a strong directional source of perturbations such as the 
Sun, which drives one component (let us say, w~^) and the 
other component is being generated at a certain distance 
from the source by means of, e.g., reflection from density 
inhomogeneities. This problem is closely related to the 
problem of the imbalance, however, in most of this paper 
we considered idealized spatially homogeneous case. Dis- 
sipation of waves in turbulence (Beresnyak & Lazarian 
2008c) is also related to this problem. 

We realize, that by relaxing homogeneity we consid- 
erably broaden the physical scope of the model, to the 
point that we might be unable to draw clear conclusions 
on the nature of inhomogeneous turbulence. In this sit- 
uation we preferred to take a first step by considering a 
toy problem of weak turbulence with two wave sources 
separated by a certain distance. 

One remarkable new property that could arise from 
such formulation is that the cascading might not reach 
the dissipation scale. This could happen if, e.g., the 
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sources of the waves are too close to each other. One 
of the questions that we are asking in this appendix, is 
whether it is possible that only one wave component, e.g. 
w" is dissipated but the other is simply distorted. 

We used simplified equations of weak cascading, the 
diffusive one dimensional k-space equations for weak per- 
pendicular cascade from Lithwick & Goldreich (2003) 
which were expanded to one spacial dimension by intro- 
ducing advection in space, making them an advection- 
diffusion equations with advection happened in real 
space, while diffusion represented energy cascading in 
k-space. We used open boundaries, so that the non- 
cascaded waves were allowed to freely escape through 
them. 

By solving those equation numerically, changing the 
distance between wave sources we found that the ap- 
proximate equality of the energies of the waves at the 
smallest scale which is reached by cascading is rather ro- 
bust feature, that is produced by cascading itself. Note 
that Lithwick & Goldreich 2003 assumed that this "pin- 
ning" is due to the dissipation term. In our toy model, 
however, regardless of whether an actual dissipation took 
place in the system, the spectra were "pinned" on small 
scales. 

Another possibility that was opened by the inhomoge- 
neous formulation was to drive wave sources with arbi- 
trary power and still obtain a stationary state. In the 
homogeneous formulation the stationary state was not 
possible if the rate of energy driving e+/e~ was larger 
than 2 (Lithwick & Goldreich 2003). Our inhomoge- 
neous toy model can deal with larger imbalances since the 
waves were allowed to escape through the open bound- 
aries of the box. Fig. [14] shows spectrum for both types 
of waves in such case of strong imbalance, and incom- 
plete cascading mentioned above (energy hasn't reached 
dissipation scale). Note, that unlike DNS with periodic 
boundaries from the main body of this paper, where en- 
ergy can only be lost through dissipation, in this case 
energy could be lost due to open boundaries and the 
stationary state could be achieved without any physical 
dissipation actually taking place. 

Fig. [Ml demonstrates a feature which was not observed 
in Lithwick & Goldreich (2003), namely that imbalance 
reverses on scales an order of magnitude smaller than the 
driving scale. This unexpected feature is fairly robust in 
the case "incomplete" cascading and strong imbalance. 
This suggest that inhomogeneous imbalanced turbulence 
could be much more complicated that its homogeneous 
counterpart and further research is necessary. 
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